Transforming Environmental Justice Research Through Open Science

# Define counties of interest
counties = ['Wayne', 'Oakland', 'Macomb']

# Fetch the county boundaries for Michigan from pygris dataset
metro_counties = pygris.counties(state="MI", year=2022)

# Filter the DataFrame to only include counties in the Detroit metro area
detroit_metro = metro_counties[metro_counties['NAME'].isin(counties)]

# Combine the geometries of selected counties into a single polygon
detroit_metro = detroit_metro.dissolve(by='STATEFP')

# Convert to GeoDataFrame and ensure coordinate reference system is EPSG:4269
detroit_metro = gpd.GeoDataFrame(detroit_metro, geometry='geometry', crs='EPSG:4269')

# Obtain the total bounding box for the Detroit metro area polygon
bbox = detroit_metro.total_bounds

# Initialize list to hold TRI facility data for each county
tri_data = []

# Loop through each county to fetch TRI data separately
for county in counties:
    # Construct the API URL for the current county
    api_url = f"https://data.epa.gov/efservice/tri_facility/state_abbr/MI/county_name/{county}/JSON"
    response = requests.get(api_url)

    # Check if the response was successful
    if response.status_code == 200:
        county_data = response.json()
        # Append the current county's TRI data to the overall list
        tri_data.extend(county_data)
    else:
        print(f"Failed to fetch data for {county} County. Status code: {response.status_code}")

# Create a Pandas DataFrame from the collected TRI data
tri_df = pd.DataFrame(tri_data)

# Reproject TRI facility data to Web Mercator for map visualization
gdf_tri_form_r_bm = gdf_tri_form_r.to_crs(epsg=3857)

# Create scatter plot with graduated symbols based on air release amount
scatter = ax.scatter(gdf_tri_form_r_bm.geometry.x, gdf_tri_form_r_bm.geometry.y, 
                     s=gdf_tri_form_r_bm['LOG_AIR_RELEASE']*20,  # Size by log of release
                     c='orangered',  # Color of symbols
                     edgecolor='yellow',  # Edge color for visibility
                     alpha=0.7)  # Transparency

import xarray as xr
import rasterio
import rasterio.mask

# Define paths to SVI TIF files - different vulnerability aspects
tif_files = [
    "data/svi/svi_2020_tract_overall_wgs84.tif",
    "data/svi/svi_2020_tract_minority_wgs84.tif",
    "data/svi/svi_2020_tract_socioeconomic_wgs84.tif",
    "data/svi/svi_2020_tract_housing_wgs84.tif",
    "data/svi/svi_2020_tract_household_wgs84.tif"
]

# Process each SVI layer
data_arrays = []
for file in tif_files:
    with rasterio.open(file) as src:
        # Clip to Detroit metro area
        metro_reprojected = detroit_metro.to_crs(src.crs)
        out_image, out_transform = rasterio.mask.mask(
            src, metro_reprojected.geometry, crop=True)
        
        # Create DataArray with coordinates
        da = xr.DataArray(out_image[0],
                         coords={'y': ('y', ys[:, 0]),
                                'x': ('x', xs[0, :])},
                         dims=['y', 'x'])
        data_arrays.append(da)

# Combine DataArrays with named layers
layer_names = ['Overall', 'Minority', 'Socioeconomic', 
              'Housing', 'Household']
svi_detroit = xr.concat(data_arrays, dim='layer')
svi_detroit = svi_detroit.assign_coords(
    layer=('layer', layer_names))

from rasterio.enums import Resampling

# Select the 'Overall' layer from SVI dataset
svi_overall = svi_detroit.sel(layer='Overall')
svi_overall = svi_overall.rio.write_crs("EPSG:4326")

# Reproject SVI to match the TRI air release raster
svi_reprojected = svi_overall.rio.reproject_match(
    air_release_raster_da)

# Transform and normalize the air release data
# 1. Log transform to reduce impact of outliers
air_release_log = np.log1p(air_release_disaggregated)

# 2. Scale log values to 0-1 range to match SVI
air_release_scaled = (air_release_log - air_release_log.min()) / \
                    (air_release_log.max() - air_release_log.min())

# Multiply scaled air release with SVI to create vulnerability index
# High values = areas with both high pollution and high vulnerability
vulnerability_indicator = air_release_scaled * svi_reprojected

# Find locations with highest combined vulnerability
vulnerability_df = vulnerability_indicator.to_dataframe(
    name='index').reset_index()
top_10 = vulnerability_df.sort_values(
    'index', ascending=False).head(10)

from rasterio.enums import Resampling

# Select the 'Overall' layer from SVI dataset
svi_overall = svi_detroit.sel(layer='Overall')
svi_overall = svi_overall.rio.write_crs("EPSG:4326")

# Reproject SVI to match the TRI air release raster
svi_reprojected = svi_overall.rio.reproject_match(
    air_release_raster_da)

# Transform and normalize the air release data
# 1. Log transform to reduce impact of outliers
air_release_log = np.log1p(air_release_disaggregated)

# 2. Scale log values to 0-1 range to match SVI
air_release_scaled = (air_release_log - air_release_log.min()) / \
                    (air_release_log.max() - air_release_log.min())

# Multiply scaled air release with SVI to create vulnerability index
# High values = areas with both high pollution and high vulnerability
vulnerability_indicator = air_release_scaled * svi_reprojected

# Find locations with highest combined vulnerability
vulnerability_df = vulnerability_indicator.to_dataframe(
    name='index').reset_index()
top_10 = vulnerability_df.sort_values(
    'index', ascending=False).head(10)

from scipy.interpolate import griddata

# Access CDC PLACES data via GeoJSON API
url = "https://data.cdc.gov/resource/cwsq-ngmh.geojson"
county_filter = " OR ".join([f"countyname = '{county}'" 
                           for county in detroit_counties])
params = {
    "$where": f"stateabbr = 'MI' AND ({county_filter})",
    "$limit": 50000
}
response = requests.get(url, params=params)
gdf = gpd.read_file(response.text)

# Filter for asthma data
gdf_asthma = gdf[gdf['measure'] == 'Current asthma among adults'].copy()
gdf_asthma['data_value'] = pd.to_numeric(
    gdf_asthma['data_value'], errors='coerce')

# Prepare for IDW interpolation
X = gdf_asthma.geometry.x.values
Y = gdf_asthma.geometry.y.values
Z = gdf_asthma['data_value'].values
mask = ~np.isnan(Z)
X, Y, Z = X[mask], Y[mask], Z[mask]

# Create interpolation grid
grid_resolution = 0.025
x_min, y_min, x_max, y_max = gdf_asthma.total_bounds
grid_x = np.arange(x_min, x_max, grid_resolution)
grid_y = np.arange(y_min, y_max, grid_resolution)
grid_xx, grid_yy = np.meshgrid(grid_x, grid_y)

# Interpolate asthma prevalence surface
grid_z = griddata(np.column_stack((X, Y)), Z, 
                 (grid_xx, grid_yy), method='linear')

from pysal.explore import esda
from pysal.lib import weights
from splot.esda import moran_scatterplot
import shapely

# Create GeoDataFrame with aligned air release and asthma data
df = gpd.GeoDataFrame({
    'air_release': air_release_aligned.values.flatten(),
    'asthma': ds_clipped.asthma.values.flatten(),
    'geometry': [
        shapely.geometry.Point(x, y) 
        for x, y in zip(
            np.repeat(air_release_aligned.x.values, 
                     len(air_release_aligned.y)),
            np.tile(air_release_aligned.y.values, 
                   len(air_release_aligned.x))
        )
    ]
})
df = df.dropna()  # Remove rows with NaN values

# Create spatial weights matrix (k-nearest neighbors)
w = weights.distance.KNN.from_dataframe(df, k=8)
w.transform = 'r'  # Row-standardize weights

# Calculate Bivariate Moran's I for air release vs asthma
moran_bv = esda.Moran_BV(df['air_release'], df['asthma'], w)
print(f"Bivariate Moran's I: {moran_bv.I}")
print(f"p-value: {moran_bv.p_sim}")

# Test shows significant spatial relationship between 
# air releases and asthma prevalence (p < 0.001)

from pysal.explore import esda
from pysal.lib import weights
from splot.esda import moran_scatterplot
import shapely

# Create GeoDataFrame with aligned air release and asthma data
df = gpd.GeoDataFrame({
    'air_release': air_release_aligned.values.flatten(),
    'asthma': ds_clipped.asthma.values.flatten(),
    'geometry': [
        shapely.geometry.Point(x, y) 
        for x, y in zip(
            np.repeat(air_release_aligned.x.values, 
                     len(air_release_aligned.y)),
            np.tile(air_release_aligned.y.values, 
                   len(air_release_aligned.x))
        )
    ]
})
df = df.dropna()  # Remove rows with NaN values

# Create spatial weights matrix (k-nearest neighbors)
w = weights.distance.KNN.from_dataframe(df, k=8)
w.transform = 'r'  # Row-standardize weights

# Calculate Bivariate Moran's I for air release vs asthma
moran_bv = esda.Moran_BV(df['air_release'], df['asthma'], w)
print(f"Bivariate Moran's I: {moran_bv.I}")
print(f"p-value: {moran_bv.p_sim}")

# Test shows significant spatial relationship between 
# air releases and asthma prevalence (p < 0.001)

from pysal.explore import esda
from pysal.lib import weights
from splot.esda import moran_scatterplot
import shapely

# Create GeoDataFrame with aligned air release and asthma data
df = gpd.GeoDataFrame({
    'air_release': air_release_aligned.values.flatten(),
    'asthma': ds_clipped.asthma.values.flatten(),
    'geometry': [
        shapely.geometry.Point(x, y) 
        for x, y in zip(
            np.repeat(air_release_aligned.x.values, 
                     len(air_release_aligned.y)),
            np.tile(air_release_aligned.y.values, 
                   len(air_release_aligned.x))
        )
    ]
})
df = df.dropna()  # Remove rows with NaN values

# Create spatial weights matrix (k-nearest neighbors)
w = weights.distance.KNN.from_dataframe(df, k=8)
w.transform = 'r'  # Row-standardize weights

# Calculate Bivariate Moran's I for air release vs asthma
moran_bv = esda.Moran_BV(df['air_release'], df['asthma'], w)
print(f"Bivariate Moran's I: {moran_bv.I}")
print(f"p-value: {moran_bv.p_sim}")

# Test shows significant spatial relationship between 
# air releases and asthma prevalence (p < 0.001)